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ABSTRACT 

We present results of a study of the virial state of high redshift dark matter haloes in 
an N-body simulation. We find that the majority of collapsed, bound haloes are not 
virialized at any redshift slice in our study (z = 15 — 6) and have excess kinetic energy. 
At these redshifts, merging is still rampant and the haloes cannot strictly be treated as 
isolated systems. To assess if this excess kinetic energy arises from the environment, 
we include the surface pressure term in the virial equation explicitly and relax the 
assumption that the density at the halo boundary is zero. Upon inclusion of the surface 
term, we find that the haloes are much closer to virialization, however, they still have 
some excess kinetic energy. We report trends of the virial ratio including the extra 
surface term with three key halo properties: spin, environment, and concentration. 
We find that haloes with closer neighbors depart more from virialization, and that 
haloes with larger spin parameters do as well. We conclude that except at the lowest 
masses [M < 10^ Mq), dark matter haloes at high redshift are not fully virialized. 
This finding has interesting implications for galaxy formation at these high redshifts, 
as the excess kinetic energy will impact the subsequent collapse of baryons and the 
formation of the first disks and/or baryonic structures. 

Key words: cosmology: dark matter galaxies: high-redshift 



1 INTRODUCTION 



For an isolated collapsed, bound dark matter halo in 
equilibrium, the scalar virial theorem, 



. In the standard A cold dark matter (LCDM) model, small 
' Gaussian perturbations in the dark matter density field at 
early times provide the seeds for the formation of struc- 
1 ture in the Universe. The dark matter gravitationally col- 
lapses, and forms bound structures that eventually relax 
into a state of virial equilibrium. The number density of 
collapsed dark matter haloes at a given mass and epoch, 
the mass function, of these relaxed haloes provides powerful 
constraints on the p arameters of the LCDM model (see e.g. 
iHaiman et al.l 1200 ll : [Cunha & Evrard 2 01Cl'). Various ana- 
lytica l methods (e.g. [Press fc Schechteij 1 19741 : 1 Bond et abl 
Il99ll : ISheth fc TormenI 1 19991 ') also predict the halo mass 
function for collapsed, bound, and virialized haloes in 
LCDM. However, simulations at high redshift {z > 1) have 
found that the majority of collapsed, bound haloes are 
not in virial equilibriu m ( Jang-Condell fc Hernguisl 200ll : 
iHetznecker fc BurkertI l2006l : iDavis fc NataraianI l201(t ). 
Thus, there seems to be a mis-match with simulations; they 
find mass functions (which assume the haloes are virialized) 
that match the analytic predictions, and yet the detailed 
structure of the haloes shows that they are not virialized. 
In this paper, we explore in detail the virialization state of 
dark matter haloes at high redshift in order to understand 
this discrepancy. 



2K + U ^0, 



(1) 



provides a simple relationship between the halo's total ki- 
netic (K) and potential (U) energies. In LCDM, dark mat- 
ter haloes are expected to reach virial equilibrium rapidly 
upon collapse when they detach from the Hubble flow. The 
timescale for virialization is of the order of the dynamical 
time, which for a dark matter halo may be estimated as 
tr ~ Ri7s/vciTc, where Rns is the virial radius and Vdic the 
circular velocity, iicirc = \/GM/Ri7s. For a 10 M© halo at 
z = & this is roughly 1 x 10* yrs, or one percent of the Hub- 
ble time at that redshift. Therefore, despite rapid merging 
activity these haloes have had sufficient time to reach virial 
equilibrium, but do not appear to do so in the simulations. 
Here we explore the energy budget of these haloes, to 
determine why simulated haloes at high redshift are ap- 
parently out of virial equilibrium. There are several pos- 
sibilities which may explain this finding. In this paper, we 
probe this issue by relaxing two assumptions typically made 
when applying the virial theorem. First, we include the non- 
negligible contributions of the environment to the halo's 
gravitational potential and secondly, we do not truncate 
the density profile of the halo at the virial radius. While 
these two assumptions are valid for isolated haloes typical 
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of the local Universe, at high redshift when the Universe 
was denser, these assumptions are incorrect. In addition to 
these two assumptions, it is possible that systematic errors 
from the halo finding algorithm may bias measurements of 
halo virialization. Finally, we explore the possibility that 
the departure from virialization is correlated with key halo 
properties. The organization of this paper is as follows: we 
derive the virial equation from the momentum equation in 
section [5] explicitly retaining the surface terms that are usu- 
ally neglected. In section[3]present simulation results and the 
correlations between halo energetics and key halo structural 
properties, namely, halo spin, concentration, and local envi- 
ronment. We conclude in section |4] with a discussion of our 
findings and their implications. 



2 THE VIRIAL EQUATION 

Here we derive the virial theorem for dark matter haloes in- 
cluding boundary effects and an external gravitational po- 
tential. An isolated system of coUisionless particles in equi- 
librium satisfies the scalar virial theorem (see Equation [T]) . 
This result is derived under the assumption that the mass 
density approaches zero at large radii for a finite mass dis- 
tribution. In this work, our systems of interest are high red- 
shift dark matter haloes, which form within a cosmic web 
of sheets, voids, and filaments. We therefore cannot assume 
that the mass density approaches zero at the boundary of 
these systems. We also consider the potential due to nearby 
particles exterior to the haloeslj 

The virial theorem is obtained from the conservation of 
momentum for a dark matter halo modeled as a coUisionless 
fluid. Following iBinnev fc Tremaind l|l987h . the correspond- 
ing fluid equation is multiplied by Xk and integrated over the 
volume V, with the caveat that the surface terms resulting 
from integration by parts do not vanish. In this case, the 
tensor virial theorem corresponding to a coUisionless fluid is 
given by. 



d% 
di2 



= 2K,J+W^J- 



pXiVjVknkdb — -—- 
2dt 



pXiXjVkTikdS, (2) 



where Uj = J pXiXjd^x is the moment of inertia tensor, 
Kij = ^ J pviVjd^x is the kinetic energy tensor, Wij = 
— j pXid^/dxjd?x is the potential energy tensor. We de- 
note the mass density, position, average velocity, and grav- 
itational potential by p, x, v, and $ respectively. The vec- 
tor n is the outward-pointing unit normal to the surface S 
that encloses the volume V. The third term on the right of 
equation ((2| corresponds to the kinetic stresses on the halo 
boundary. The last term is the time d erivative of the mo- 
ment of inertia flux through the surface (|Ballesteros-Paredesl 
120061 ). The scalar virial theorem is obtained by taking the 
trace of equation ([2]) : 



d^J 



px ■ V V ■ dS — -— (P px^ V ■ dS. (3) 



FoUowing iBiiinev fc Tremaind (| 19871 ) and 

I Ballesteros- Pareded (|200q ). the potential energy term 



^ In what follows, we use the Einstein summation convention 
where repeated indices correspond to summations. 



may be broken up into contributions from particles inside 
and outside of the halo, $ = "I>int + ^cxt , so that 



W 



1 



P 'I'intdV 



p X- dV. 

ox 



(4) 



Note that the flrst term is the gravitational potential energy 
of the halo, which we denote as U from here on. In the case of 
a spherical volume, the last term in equation Q is B?M /2, 
where R is the radius and M is the mass enclosed within the 
volume. Assuming that this term is negligible for a steady- 
state system, we obtain 



= 2K + U + U,^t- E. 

t/oxt = 



- 9$cxt ,,, 
P X ■ — ;^^— dV, 



Es = 



px ■ V V 



dx 



■dS, 



(5) 
(6) 

(7) 



for a halo in virial equilibrium. Ucxt and Es represent cor- 
rections to the standard virial relation and can be explicitly 
calculated for simulated dark matter haloes. 



3 SIMULATION RESULTS 



Using GADGET-2 (|Springel|[2005l '). we ran a dark matter 
only simulation with the WMAP5 cosmological parameters 
(|nM, nA,nb,/i,n,a8} = {0 .258, 0.742, 0.044, 0.719, 0.963, 
0.7961. iDunklev et all 120091 ') from 2 « 100 to 2 = 6 (see 
iDavis fc NataraiarJ 1 20101 ) for full details). Our dark matter 
particle mass is m = 1.0 x lO^Mo/h, which sets a comov- 
ing box size at 2.46 Mpc/h, with 512^ dark matter particles. 
We use the HOP algorithm to iden tify collapsed dark mat- 
ter haloes (|Eisenstein fc Hutill998l ). HOP first calculates a 
density for each particle by smoothing over its nearest 64 
neighbors using a cubic spline kernel. It then groups parti- 
cles with their densest neighbor, and density thresholds are 
used to ensure that haloes are not being over-counted as 
sub-haloes within a larger halo. We choose density thresh- 
ol ds to match the high redshift mass function described 
in iReed et al.l l|2007r ). Lastly, we remove unbound particles 
from each halo. 

For each halo, we calculate the total potential energy, 
U, using a direct summation over the particles assigned to 
the halo: 



U 



•sr-^ \-^ Grriimj 



E 



i — 1 j^i + 1 



(8) 



where rij is the separation between particles i and j, G is, 
Newton's gravitational constant, and m the particle mass. 
The total kinetic energy, T is found by summing over the 
particle's individual kinetic energies: 



N 
1 ^"^ _,2 



(9) 



In measuring T and [/, we must ensure that we 
have enough particles in our halo sample to accurately 
measure the ener gi es. W e addressed this question in 
IPavis fc NataraiarJ ([20101), and summarize the relevant 
points here. We re-ran our simulation with one eighth as 
many particles (256'' total) and with eight times as many 
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Figure 2. Scaled histogram of /9 evaluated at four different radii. 
We find that the cores of our haloes are in virial equilibrium, 
while the overall halo is not. 
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Figure 1. Fractional difference of T (top) and U (bottom) as 
a function of the number of particles in the low resolution halo. 
The error bars represent the la deviation about the mean. We 
find that at least 300 particles are required to accurately measure 
the kinetic and potential energies in our sample. 



particles (1024'' total). This allows us to compare directly 
individual halos with different resolutions. We used the halo 
centre of mass to cross-matched the halo catalogues. We 
show in Figure[T]the change in T and U for the 512'' and 256^ 
runs as a function of the number of particles in the low res- 
olution halo. We find that at least 300 particles are required 
to accurately measure the kinetic and potential energies, and 
restrict our halo sample to halos with at least this many par- 
ticles, corresponding to a halo mass of M = 3 x 10^ Mo/h. 
Having measured the kinetic and potential energies of 
our haloes, we can define the virial ratio. 



I3^2T/U + 1, 



(10) 



which for fully virialized haloes should be zero. However, 
as our measurements of T and U are instantaneous quan- 
tities, and not time averaged, we expect a range of values 
for /3, with a mean value of zero. There are various cuts on 
P used in t h e lite rature to select st rictly virial ized haloes: 
IShaw et al.l (|2006l) use P > -0 .2, iBett et all (J2007l ) use 



< 0.5, and iNeto et al.l hOOJl ) use /3 < 0.35. As noted 
Davis fc NataraianI ||201(J ) we find that our haloes do not 



1/3 
in 

have (/3) « 0, but are offset to high values of kinetic energy 
such that (P) < for all redshifts, and haloes are further 
from virialization at higher redshifts. One possibility is that 
the extra terms in the virial equation (C/oxt and Ea) are 
not negligible at higher redshifts, where large amounts of 
infalling material contribute to the terms. As the Universe 



expands, we expect a smaller contribution from these ex- 
tra terms. Our finding s follow the general trend reported in 
iHetznecker fc Burkerd pOOd ) at lower redshift, though their 
fitting function cannot be extended to our redshift range. 

The virialization process is expected to happen from the 
inside out. Thus, it should be possible to identify a virialized 
core. To do this, we calculated T and U only for particles 
inside three smaller radii: 0.757?i78,0.57?i78, and Q.257?i78. 
For the choice of boundary at 0.25i?i78 , we do find virialized 
cores, as expected. We show in Figure [2] histograms of /? for 
these three inner radii as well as the histogram of /3 for the 
entire halo. We find that for smaller radii, the mean value of 
P shifts towards 0, as expected for a virialized object. Thus 
we conclude that the halo cores of our sample are virialized, 
while the entire halo is not. 

The departure from virial equilibrium motivated our in- 
terest in the extra terms in the virial equation. We calculate 
the contribution to the potential energy from the external 
mass distribution, f/cxt (equation [S}. Converting the volume 
integral into a sum over particles, we find 



Ucxt — 



E 



Tfli Xi • dij 



(11) 



where A'^ is the total number of particles in the halo (see 
equation llip . To calculate the gravitational acceleration a 
acting on particle i we sum over all particles, j, within 
107?i78 which are not part of the halo itself: 



JVoxt 



i 



GrajVi- 



(12) 



However, after calculating this term, we find that it is neg- 
ligible compared to U (of order 1%). Thus, even at higher 
redshifts where it may be expected that external matter will 
affect the total gravitational term in the virial equation, we 
find that it provides only a small contribution to the total 
energy. 

We next exam ine the term in equation [71 Eg. Following 
IShaw et al.l l|2006l ) , we select halo particles between 0.87?i78 
and 7?i78, where 7?i78 is the radius which encloses a mean 
overdensity of 178pcrit. This defines the surface over which 
we perform the summation. Then, approximating the sur- 
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face integral as a sumiiiation over particles in this shell, we 
calculate Ea as: 



Es 



Rn 



Ar 



■E" 



(13) 



where Rm is the mean radius of the Ns particles in the shell 
of thickness Ar = 0.27?i78- Each particle, i, has a radial ve- 
locity Vi,r = fi-Vi/\fi\. We note that our Ea term is similar to 
that of lShaw et alj (120061 ). but not identical to their surface 
pressure term. Here we included only the contribution from 
the velocity component normal to the surface. This may be 
important if the velocities of the particles falling into the 
halo are on predominantly radial orbits, such as falling into 
the halo along a filament, and would then have a larger sur- 
face term than when simply calculating the total pressure 
(which is proportional to |t/p/3). We find that at least 30 
particles are required in the shell to correctly resolve and 
calculate Es; haloes with fewer particles have strongly bi- 
ased (low) values of Es- We therefore only include haloes 
with more than 30 particles in this shell for the remainder 
of this work. 

We point out that our halos are not spherical - there is 
no guarantee that the density cuts used by the HOP algo- 
rithm yields spherical halos. However, in our calculation of 
Ea, we have chosen a spherical surface over which to measure 
this term. Thus, we wish to test whether there is any bias 
in the value of Ea arising from the real shapes of haloes. 
To measure halo shapes, we use the same method as in 
iDavis &: NataraianI l|20ld ). First, we calculate the normal- 
ized moment of inertia tensor. The ratios of the eigenvalues 
(a > fe > c) of this tensor can be used to define a halo's 
sphericity (s = c/a) and triaxiality (t = (a^ — (?)/{a? — b^)). 
We found, however, no trend between shape and Es or /3. 
Thus, we conclude that our choice of using a spherical shell 
rather than following the outer boundary of the HOP halo 
does not bias our results. 

Having found that [/ext is negligible, we now include 
only the surface pressure term, Ea, in our calculation of a 
corrected virial ratio, /?', where: 



P' = {2T^Ea)/U + l, 



(14) 



and show in the top panel of Figure [3l a histogram of 
/3 and /3' at z = 6. Without including the Es term, we 
find (/3) = -0.295 ± 0.0010, and after inclusion, we find 
(/3') = -0.076 ± 0.0013. The standard deviation of the cor- 
rected distribution increases slightly, from cr^ = 0.097 to 
a pi = 0.126. With the inclusion of Es, we find that the 
virial theorem is nearly satisfied by our haloes a t 2: = 6: ap- 
proxim ately 84% lie within the cut applied by IShaw et al.l 
i|2006h of /3' > -0.2, whereas only 15% of the haloes have 
/3 > -0.2. 

At higher redshifts, however, {/?') becomes increasingly 
negative, implying that our haloes are farther from virial- 
ization at higher redshift even after the correction arising 
from the surface term is taken into account. This is to be 
expected, as these haloes have had less time to fully relax. 
We show in the bottom panel of Figure [3] (/3) and {/?') as a 
function of redshift. We include as the so lid curve the fitting 
function of lHetznecker fc BurkertI ((20030 which they found 
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Figure 3. Top panel: histograms of the virial ratio, /9, and the 
corrected virial ratio, /3', at 2: = 6. The mean and standard de- 
viation of each distribution is also shown. Bottom panel: mean 
virial ratio, (/3) (diamonds), and mean corrected virial ratio, (/3') 
(asterisks), as a function of red shift. The solid curve is the fi tted 
polynomial for (/3) provided bv lHetznecker fc BurkertI 1120061 ) for 
redshifts 2 < 3. The error bars represent the la deviation about 
the mean. 



for redshifts less than 3. We find that the correction term, 
Ea does not fully correct the virial ratio at any redshift, and 
the in fact correction is less effective at higher redshifts as 
seen in Figure [3] 



3.1 Correlation of /?' with halo properties 

After including the extra term, Es from the virial equation, 
we have found that the haloes still do not have a mean /3' of 
zero as expected. We compare /?' to other halo properties, 
in an attempt to find the source of the offset. 

The effect of the local environment is first examined, 
as it addresses the assu mption of isolation in the v irial the- 
orem. In a recent paper iDavis fc NataraianI (|201(J ) studied 
correlations between the local environment and halo struc- 
tural properties. We use one of their metrics of the local 



Note that there is a typo in their fitting function as reported 



in their paper. It should read 77(2) 
where r] = 2T/\U\. 



-3.3x10^2(2-2.7)2-1-1.3, 
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environment, the distance to the 3'^'* nearest neighbor (D3), 
and compare it with the virial ratios, /3 and /?'. We expect 
that haloes in denser environments and with close neigh- 
bors will have larger surface terms, as these haloes are least 
likely to be in isolation as is assumed by the standard virial 
theorem. Thus we expect a trend with P and environment, 
and after accounting for the surface terms, any remaining 
correlation could explain the increased dispersion seen in /3' 
compared to p. 

We find that for extremely unvirialized haloes such that 
/3 < —0.45, the mean value of D3 is 5.25 kpc with a standard 
deviation of cr = 1.92 kpc. We note here that our box size 
at this epoch is ~ 580 kpc. However, for the full sample of 
haloes, the mean distance is 6.19 kpc with a = 2.94 kpc. 
Thus, the extremely unvirialized haloes are slightly system- 
atically closer to other haloes. In Figure |3] we show the cor- 
relation between /?' and D3 at z = 6; the solid curve shows 
the mean value of D3 binned by /3' with error bars depict- 
ing the \a error of the mean, and the dashed curve shows a 
linear fit to the mean trend with a slope m — 2.51 and y- 
intercept b = 6.39. The fit has a reduced x^ = 1-78. There is 
only a slight trend for haloes with small values of P' to have 
closer neighbors than halos with large values of /3'. However, 
there is large scatter in the values of D3 in a given range of 
/3', implying that environment is not the only factor in the 
dispersi on of /?'. 

Hetznecker fc BurkertI l|2006l ) report a correlation be- 
tween the spin parameter, 



A = 



J\T + U\ 



1/2 



and their definition of the virial coefficient, rj — 2T/\U\. 
They found that for redshifts between — 3, A is propor- 
tional to J?*, and that the relationship varied with redshift. 
As they note, we expect a relationship between A and /3', 
with scatter due to the angular momentum, J, in the defini- 
tion of A. When we use their definition of r/, we find similar 
results, but this definition does not account for the surface 
term of equation 1131 We show in Figure |4] the correlation 
between A and P' ai z — 6. We report a similar trend to 
iHetznecker fc BurkertI (|2006l ). For comparison, we fit our 
data to the same relationship: A = a + ^{P' — l)**, and found 
best fit va lues of a = 0.04 and 7 = 2. 6 x 10~^, similar to 
the values IHetznecker fc BurkertI l|2006r ) found of a = 0.036 
and 7 = 2.4 x 10"'^ for all haloes in their simulation between 

< z < 3^ 

In iDavis fc NataraianI (|201G| ') we fit NFW profiles to 
our halo sample at z = 6. Using the scale radius from the 
fit, we define the halo concentration, Cns = Rns/vs- In 
this work, we look for a correlation between halo concen- 
tration and P' . Analytic models of dark matter haloes with 
galaxy and cluster scale masses at low redshift and which 
use the NFW density profile find that haloes with smaller 
concentrations a r e furt her from virialization (/3 < —0.2) 
iLokas fc Mamon 1I2OOH). F or these massive haloes, simula- 



Shapiro et al.l ||20 04') match the analytic trends 



Lokas fc MamonI (|2001|). We find, however, that 



tions run by 

predicted by 

our halo sample d oes n ot follow the same trend reported in 

iLokas fc MamonI ||200H ). Our simulations show that haloes 

which are virialized (/3 > 0) have smaller concentrations 

than haloes with P < —0.2. 

Another method of finding relaxed haloes is to use the 



offset between the center of mass and the cent er of the po- 
tential well (see e.g. lD'Onghia fc Navarroll2007h . We use the 
location of the densest particle in the halo as a proxy for 
the center of the halo's potential well, and compare this off- 
set. As = (a;cM — a;dcn)/i?i78, to P' . We find a small trend 
for haloes with small values of /?' to have larger offsets, as 
shown in Figure |4] for the haloes at z = 6. However, we note 
that there does not seem to be a systematic shift in /3' for 
haloes with As < 0.1. Thus, we conclude that using As to 
find virialized haloes is not effective at these redshifts. 



4 DISCUSSION 

In this paper, we explore the question of whether or not 
high redshift dark matter haloes are virialized. Including 
the surface term in the virial equation, 82% of our haloes at 
z — 6 are virialized (/?' > —0.2) whereas without the term, 
only 15% of them are. However, the distribution of virial 
ratios is still not centered around 0. A mean of would be 
expected if the spread of /3' is due simply to the fact that 
our measurements of the kinetic and potential energies are 
not time-averaged. Therefore, we conclude that simulated 
dark matter haloes at high redshifts are not virialized. On 
average, they have too much kinetic energy. One explanation 
for this excess kinetic energy is due to the approximations 
that we have made in calculating the terms in the virial 
theorem. A key issue is the choice of the shell in calculating 
Ea, which is arbitrary. It needs to be large enough so that we 
can sufficiently resolve the flux through the shell. However, 
with better particle resolution, the shell could be narrowed 
for a better approximation to the surface integral. We have 
also flxed the shell as spherical, whereas the haloes are not 
truly spherical. Also, because of the time resolution of our 
snapshots, we are not able to test if the moment of inertia 
fiux (the last term in equation [Sl is truly negligible. 

Our findings could be a consequence of the deployed 
halo finding algorithms at high redshifts. Since virialization 
happens from the inside out, it is possible to find a virialized 
core by looking at regions interior to the standard definition 
of the virial radius. This would be necessary if one seeks to 
find fully virialized sets of particles at high redshift. Alter- 
ing the halo finder would however, alter the mass function 
derived from our simulations, which would put it out of line 
with mass functions fro m other simulations at similar red- 
shifts (|Reed et al.ll2007l ') and from analytic predictions. 

Finally, our finding has important implications for 
galaxy formation. Haloes with excess kinetic energy will 
mostly likely also have increased turbulence when baryonic 
gas falls into the central regions of the halo. The ensuing tur- 
bulent viscosity could be an effective method for transport- 
ing angular momentum in high redshift dark matter haloes. 
This might alter the cooling and collapse of baryons in these 
high redshift haloes than expected, as suggested in recent 
results reported by Grief et al. (2011). 
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Figure 4. Correlatio ns of distance to third nearest neighbor (D3, top left panel), spin parameter ( top right panel, plotted versus (/3' — 1)^ 
for comparison with iHetznecker fc BurkertI 1 120061 )). concentration (Ciyg, bottom left panel), and offset parameter (As, bottom right) 
with the corrected virial ratio, /?', for the haloes in our sample at z = 6. The contours enclose 90%, 70%, 50%, 30% and 10% of the haloes, 
and the asterisks represent mean values for each halo property. We find that haloes with closer neighbors, large spin parameters, and 
extremely large offsets are more likely to be further from virialization. We do not find the general trend reported in lShapiro et al.l ||200J) 
for concentration versus 0' . 
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